#include "rectMesh.hpp"

namespace hms {

/**
 * @brief Helper function to create ghost layer coordinates
 * 
 * @param oldCoords 
 * @return ArrayXs Inner coordinates are preserved, ghost layer coordinates
 * are added on either end
 * 
 * Repeats the last cell width for all ghost cell layers.
 * Mirroring version is implemented, but commented out.
 */
ArrayXs fillWithGhosts( const Ref<const ArrayXs>& oldCoords ){
	Index size { oldCoords.size() };
	/* reserve space for ghost cell coordinates on either Side */
	ArrayXs coords ( size + 2 );
	/* put old coordinates in the middle */
	coords.segment( 1, size ) = oldCoords;
	/* repeat last distance for ghost cell layers */
	coords(0)      = 2*oldCoords(     0) - oldCoords(     1);
	coords(size+1) = 2*oldCoords(size-1) - oldCoords(size-2);
	return coords;
}


RectMesh::RectMesh(
	ArrayXs&& xCoords,
	ArrayXs&& yCoords,
	bool ghostsIncluded
){
	if (ghostsIncluded){
		m_xCoords = std::move(xCoords);
		m_yCoords = std::move(yCoords);
	} else {
		m_xCoords = fillWithGhosts( xCoords );
		m_yCoords = fillWithGhosts( yCoords );
	}
	storeDiffs();

	m_edges.storeIndices();
	m_cells.storeData();
}


/**
 * @brief Returns the x-coordinates belonging to mesh vertices, but also those
 * used to calculate ghost cell properties.
 * @see RectMesh::xCoordsInner()
 * 
 * @return const ArrayXs& 
 */
auto RectMesh::xCoords() const -> const ArrayXs& { return m_xCoords; }
/**
 * @brief @see RectMesh::xCoords()
 * 
 * @return const ArrayXs& 
 */
auto RectMesh::yCoords() const -> const ArrayXs& { return m_yCoords; }


/**
 * @brief Returns the x-coordinates belonging to mesh vertices.
 * @see RectMesh::xCoords()
 * 
 * @return ArrayXs::ConstBlockXpr 
 */
auto RectMesh::xCoordsInner() const -> Eigen::VectorBlock<const ArrayXs> {
	return m_xCoords.segment( nGhostLayers(), nCellsX() + 1 );
}

/**
 * @see RectMesh::xCoordsInner()
 * 
 * @return Eigen::VectorBlock<const ArrayXs> 
 */
auto RectMesh::yCoordsInner() const -> Eigen::VectorBlock<const ArrayXs> {
	return m_yCoords.segment( nGhostLayers(), nCellsY() + 1 );
}

/** @brief Returns the width (length in x) of all cells including ghost cells
 * as an array. */
auto RectMesh::dx() const -> const ArrayXs& {
	return m_dx;
}
/** @brief Returns the height (length in y) of all cells including ghost cells
 * as an array. */
auto RectMesh::dy() const -> const ArrayXs& {
	return m_dy;
}

/** @brief Returns the width of all domain cells as an array 
 * @see RectMesh::dx()
 */
auto RectMesh::dxInner() const -> Eigen::VectorBlock<const ArrayXs> {
	return m_dx.segment( nGhostLayers(), nCellsX() );
}

/** @brief Returns the height of all domain cells as an array
 * @see RectMesh::dy()
*/
auto RectMesh::dyInner() const -> Eigen::VectorBlock<const ArrayXs> {
	return m_dy.segment( nGhostLayers(), nCellsY() );
}


auto RectMesh::nCellsXImp() const -> Index {
	return m_dx.size() - 2*nGhostLayers();
}
auto RectMesh::nCellsYImp() const -> Index {
	return m_dy.size() - 2*nGhostLayers();
}

auto RectMesh::sizeImp() const -> Array2i {
	return { this->nCellsX(), this->nCellsY() };
}

auto RectMesh::vertexImp(Index i) const -> Vector2s {
	return {
		m_xCoords[ i % (nCellsX()+1) + nGhostLayers() ],
		m_yCoords[ i / (nCellsX()+1) + nGhostLayers() ]
	};
}

auto RectMesh::domainAreaImp() const -> scalar {
	/* difference between coordinates of first and last vertex are the length
	 * and width of the mesh */
	return ( vertex( nVertices()-1 ) - vertex(0) ).prod();
}

void RectMesh::storeDiffs(){
	m_dx =
		m_xCoords.tail(m_xCoords.size() - 1) -
		m_xCoords.head(m_xCoords.size() - 1);
	m_dy =
		m_yCoords.tail(m_yCoords.size() - 1) -
		m_yCoords.head(m_yCoords.size() - 1);
}

}